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Abstract 

A general framework for goal-oriented a posteriori error estimation for 
finite volume methods is presented. The framework does not rely on re- 
casting finite volume methods as special cases of finite element methods, 
but instead directly determines error estimators from the discretized fi- 
nite volume equations. Thus, the framework can be applied to arbitrary 
finite volume methods. It also provides the proper functional settings to 
address well-posedness issues for the primal and adjoint problems. Nu- 
merical results are presented to illustrate the validity and effectiveness of 
the a posteriori error estimates and their applicability to adaptive mesh 
refinement. 

1 Introduction 

Finite volume methods have become increasingly popular due to their intrinsic 
conservative properties and their capability in dealing with complex domains; 
see, e.g. [11, 15]. Hence, a posteriori error estimates for finite volume methods 
are important as they aid in error control and improve the overall accuracy 
of numerical simulations. A posteriori error estimates also play a key role in 
the implementation of adaptive mesh refinement methods. In fact, the main 
motivation of the current work is to find simple and robust a posteriori error 
estimators to guide adaptive mesh refinements for finite volume methods in 
regional climate modeling [7, 16]. 

The literature on a posteriori error analysis for finite volume methods is slim 
compared to that for finite element methods. A large part of the available works 
aim to derive a posteriori error bounds for approximate solutions in certain 
global energy norms; see, e.g. [1, 2, 3, 5, 12, 13, 18, 19]. 

We are particularly interested in another type of a posteriori error estimates 
that are goal oriented, i.e., estimates of errors in certain quantities of interest. A 
goal-oriented error estimate is potentially very useful in assisting in error control. 
However, the literature on goal-oriented a posteriori error estimates for finite 
volume methods is even scarcer, probably due to the fact that finite volume 
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methods do not naturally fit into variational frameworks. Insightful efforts 
have been made to address this challenge by exploiting the equivalence between 
certain finite volume methods and numerical schemes in variational forms such 
as the finite element methods or discontinuous Galerkin methods. For example, 
in [4], a goal-oriented a posteriori error estimate is presented for a special type 
of finite volume method that is equivalent to a Pet rov- Galerkin variant of the 
discontinuous Galerkin method. In [9], an a posteriori error analysis is presented 
for ccll-ccntcrcd finite volume methods for the convection-diffusion problem by 
utilizing the equivalence between the finite volume methods and the lowest- 
order Raviart-Thomas mixed finite element method with a special quadrature. 
However, the applicability of this approach is limited for two reasons. First, 
on many grids on which finite volume schemes are constructed, there are no 
quadrature rules known, and thus implementing finite clement or discontinuous 
Galerkin schemes on these grids is impractical. An example of such grids is 
hexagonal Voronoi grid ([?]). The other reason is that finite volume methods 
for real-world problems are often sophisticated in themselves, and it is often not 
clear, to say the least, how to establish a connection to schemes in variational 
forms. In this regard we again refer to [15]. 

In this work, we aim to derive a general functional analytic framework for a 
posteriori error estimation for arbitrary finite volume methods. The idea is to 
derive a posteriori error estimators at the partial differential equation level in an 
appropriate functional setting. This approach does not require the differential 
equations or the numerical schemes to be recast in variational form nor do they 
rely on connecting an finite volume method with a finite element method. The 
approximate solutions produced by finite volume methods are simply taken 
as inputs to the a posteriori error estimator. Because the a posteriori error 
estimation is independent of the exact form of the finite volume method, it can 
be applied to arbitrary finite volume methods. 

In Galerkin finite element or discontinuous Galerkin methods, the difference 
u — Uh between the exact solution u and the approximate solution Uh is orthog- 
onal to the test function space Vh- For this reason, a posteriori error estimation 
for these methods requires that the adjoint solution </> be sought in a space Vh> 
larger than Vh- This restriction does not apply in our approach due to the 
very fact that finite volume schemes do not naturally fit into variational forms, 
though for the sake of accuracy, it may be advantageous to seek the adjoint 
solution with a higher order scheme. This point will be made clear in the next 
section. 

It has been pointed out that the well-posedness issue for the adjoint equation 
remains challenging and open in many cases; see, e.g., [4]. A byproduct of our 
approach is that, because the adjoint problem is naturally posed in an appro- 
priate functional setting, its well posedness can be dealt with by the abundant 
analytical tools of standard partial differential equation theories; again, see [4]. 

The rest of the paper is organized as follows. At the beginning of the next 
section, we present our approach for a posteriori error estimation for finite vol- 
ume methods in a general functional analytical framework. It is followed by a 
numerical example demonstrating the error estimates. In Section 3, we present 
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an application of the a posteriori error estimation to adaptive mesh refinement. 
We conclude with some remarks in Section 4. 



2 A posteriori error estimates for finite volume 
methods 

2.1 Abstract framework 

Let H denote a Hilbert space endowed with the norm || • || and inner product 
(•,•). Let L denote an unbounded operator in H with domain D(L) dense in 
H. The primal problem we deal with is succinctly formulated in this functional 
setting as: 

for each f G H, find u G D(L) such that 

Lu = f. (1) 

We consider time-dependent problems so that the operator L usually takes the 
form 

where A represents a linear differential spatial operator that is usually also 
unbounded in a respective function space. In this work, we assume that the 
primal problem (1) is well-posed, i.e., it possesses a unique solution. 

The quantity of interest is given as a possibly nonlinear functional Q(u) of 
u. Let denote an approximate solution of the primal problem (1). Then, 
the error in the quantity of interest can be written as 

Q(u)-Q(u*) = (u-u#,<f>) (2) 

for some kernel function 4> G H. We note that the kernel function <f> may depend 
on u; indeed this is the case when Q is a nonlinear function of u; see Section 3. 

The foregoing assumption that the domain D(L) is dense in H is crucial as 
it allows us to rigorously define the adjoint operator L* of L and its domain 
D(L*). Indeed, according to [17], a function u G H belongs to D(L*) if and only 
if the mapping D(L) : u — > (Lu, u) is a linear bounded functional on D(L) 
for the norm of H. Because D(L) is dense in H, the Hahn-Banach theorem 
guarantees that the bounded linear functional can be extended to the whole of 
H. Then, by the Riesz representation theorem, there exists a unique element, 
denoted as L*u, of H such that 

(Lu, u) = (u, L*u) VueD(L). (3) 

Thus, L* is a linear, possibly unbounded, operator from D(L*) to H. The 
adjoint problem can be formally stated as: 

for a given function <f> G H , find u G D(L*) such that 

L*u = (f>. (4) 
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We will demonstrate, through examples, how to compute adjoint operators. 
Generally speaking, for time-dependent problems, the operator L* usually has 
the form of 



Primal time-dependent problems are usually initial-boundary value problems. 
Consequently, most adjoint problems are /ina^-boundary value problem because 
the values of the unknown function are imposed at the final time t — T. Apply- 
ing the change of variable r = T — t transforms the final- value problem into an 
initial-value problem. In this way, many analytic techniques can be employed to 
establish the well-posedness of the adjoint problem. We will demonstrate this 
point with specific examples. 

Now that all the necessary functional settings have been introduced, we shall 
derive an a posteriori error estimate for the forward problem (1), regardless 
of the finite volume methods actually used to solve the primal and adjoint 
problems. By (2)-(4), we infer that 

Q{u) -Q{u*) =(u-u*, <P) 

={u-u*, L*u) 
={L(u - tt # ), u) 
={f-Lu*,u). 

Therefore, we have 

Q{u) - Q{u*) = (/ - Lu*, u). (5) 

Note that the true solution u is not required for computing the error and that, 
if the adjoint solution u is exact, then the error estimate is actually exact as 
well. These are the key advantages of the a posteriori estimate (5). The cost 
incurred is, of course, the need to solve the adjoint problem (4). 

So far we have not discussed the exact form of the finite volume methods. 
The approximate solutions produced by these schemes are taken as input to (5). 
Hence, the formula, in principle, applies to arbitrary finite volume methods. We 
should also note that the approximate solution to the primal problem is usually 
given in a discrete form and thus the term Lu# in (5) is not well defined at this 
point. We will explore a walk-around to this issue when we deal with specific 
examples in Section 2.2 and 3. 

2.2 Numerical demonstration: one-dimensional scalar equa- 
tion 

The primary goal of this section is, by the means of a simple example, to demon- 
strate the implementation of the abstract framework laid out in the previous 
section. We will also demonstrate the validity of the a posteriori error estimates 
by comparing the estimated errors with the true errors. Finally, we will explore 
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the impact of the numerical errors in the adjoint solution on a posteriori error 
estimates, and what measures can be taken to ensure accuracy. 
We consider the one-dimensional linear transport equation 

u t + au x = f, < x < 1, t > 0, (6) 

u(0,t)=g(t), t>0, (7) 

u(x, 0) = u (x), < x < 1. (8) 

We assume that the coefficient a is positive and constant, in which case the 
initial and boundary value problem (6)-(8) is well-posed, provided that the 
problem is cast in an appropriate functional setting. The well-posedness of this 
problem will not be discussed here because it is a classical example in partial 
differential equation theory; see, among many texts, [10]. Nevertheless, we shall 
specify the proper functional settings for the problem so that its adjoint problem 
can be defined and an a posteriori error estimate for its solution can be derived. 

We let H = L 2 ((0, 1) x (0, T)) and let L be the linear operator associated 
with (6)- (8). We define the domain D(L) as 

V(L) = {ue H \u t + au x e H, u(0, t) = 0, u(x, 0) = 0} . 

It can shown that D(L) is dense in H. For each u <E T>(L), the operator L is 
defined as 

Lu = u t + au x . 

We now define the adjoint operator L* of L and its domain V(L*). We recall 
that the domain V(L*) is defined as a space of functions u for which (Lu,u) is 
a continuous functional for u e T>(L) with respect to the norm of H , that is, 

u e V{L*) -i==^> u — > (Lu, u) is continuous in the H norm. 
Based on this definition, we determine that V(L*) is given by 

V(L*) = {u e H \u t + au x e H, u(l, t) = 0, u(x, T) = 0} 
and, for each u G V(L*), we have 

L*u = —Ut — au x . 

For u e T>(L) and u e T>(L*), the following important relationship holds: 

(Lu, u) = (u, L*u). (9) 

Within this functional setting, the adjoint problem can be stated as follows: 
for any (f> £ H, find u e D(L*) such that 

L*u = <j). (10) 
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We now consider the well-poscdncss issue for the adjoint problem (10). The 
equations of this problem, in differential form, read 



—Ut - au x = 4> 
u(x,T) = 
u(l,t) = 



for < x < 1, < t < T, 
for < x < 1, 
for < t < T. 



(11) 
(12) 

(13) 



The conditions (12) and (13) are due to the requirement that v should be sought 
in the adjoint domain T>(L*). The system (11)— (13) needs to be solved backward 
in time because a final condition is given at t = T. To overcome this awkward- 
ness, we make the change of variable t = T — t. Then, the system (11)— (13) 
becomes 



which can be solved forwards in time. The well-posedness of this system can be 
studied in the same way as that for the primal system (6)-(8), e.g., by semigroup 
theory. In the language of that theory, the operator L* is the infinitesimal 
generator a semigroup of contractions S(t). The well-posedness result for the 
system above is stated in the following theorem. 



Theorem 2.1. If <f> £ L x (0, T; i 2 (0, 1)), then there exists an unique solution 
ueC{[0,T];L 2 (0,l)) of (14)-(16) of the form 



We do not prove the theorem here and instead refer to [14]. 

The simple example considered in this section demonstrates a key advantage 
of our approach: we derive and use an adjoint system whose well-posedness can 
be studied in the same manner as that for the primal system, for which a 
multitude of analytical tools are available. 

We now concern ourselves with linear quantities of interest that can be ex- 
pressed as 



where is a kernel function independent of the unknown u. Identifying 4> i n 
(17) with <f> in (14), the a posteriori error estimate of Q(u) — Q{u#) is then 
given by (5). However, there is an implementation issue associated with (5). In 
practice, we only have the approximate solution u# in its discrete form, e.g., 
as approximations of the function values at discrete grid points. For (5) to 
make sense, we could either define a discrete analogue L# of L or define and 
apply a mapping that maps u& from its discrete representation to a continuous 
representation so that the differential operator L can be applied. In this work, 



u{x,Q) = 0, 
u(1,t)=0, 



u T — au x = 



(14) 
(15) 
(16) 





(17) 



G 



we take a third approach. We find a way to transpose L back onto u and then 
use (10) to replace L*u by <p. This approach avoids approximating Lu # with 
discrete values which usually introduces another source of error. However, we 
note that, in general, u# does not belong to T>(L) so that we cannot invoke (9) 
for (Lu#,v). Thus, we proceed by integration by parts: 



(Lu*,v)— J (uf + auf)v 



u*v — a I vrv + {u* , 4>). 

t=0 Jx=0 



Therefore, 



Q(u)-Q(u*) = (f,v)-(u*,<f>)+ [ u v + af gv. (18) 

Jt=0 Jx=0 



2.2.1 Numerical results 

The forward problem (6)-(8) is solved with the first-order explicit upwind 
method. The exact form of the scheme is not essential to this work, but we 
shall briefly state the method for the sake of reference. (Other finite volume 
methods used in the sequel, however, will not be explicitly described.) We 
choose = x n < x\ < ■ ■ ■ < Xm = 1 and let Ki = [xi-i, Xj\ denote the i th 
control volume/cell. Obviously, {K i \^L 1 forms a partition of the interval [0, 1]. 
We also set At = T/N and let t n = nAt denote the discrete time steps. We let 
{/" denote the average of the unknown u over the cell K t at time t — t n . Then, 
the first-order explicit upwind method for (6) can be written as 

jjn+l jjn Tjn _ jjn 

U * , * +a Ui s l - X =F? for Ki<M. 
At Si 1 ~ ~ 

In the above, Si = \Ki\ and Uo is the average of the unknown u on an artificial 
cell K = [x-i, x ] with a;_i = —xi. The boundary condition (7) can be 
imposed as 

To demonstrate the validity of the a posteriori error estimate (18) for the 
problem (6)-(8), we consider an example with a sine wave solution: 

u t + u x = 0, < x < 1, < t < 0.5, (19) 
u(x,0) = sin(27rx), (20) 
u(0,i) = -sin(2a7rf). (21) 

This problem has an analytic solution 

u(x, t) = sin(27r(a: - at)), (22) 

which allows us to study the performance of the a posteriori error estimate (18) 
by comparing the estimated errors to the true errors. 
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For the kernel function in (17), we first consider the simple case 

0=1. (23) 

This choice of kernel function emphasizes the accuracy of the solution every- 
where in the spatial-temporal region [0, 1] x [0, T). With this choice of kernel 
function, we compute the solution on an array of uniform grids, from a coarse 
one with M = 20 * 2° cells to a fine one with M = 20 * 2 4 cells. Because the 
solution of (19)-(21) is just the sine wave given by (22), we can compute the 
true error in the quantity of interest Q(u) — J Q udxdt. We also compute the 
solution of the adjoint equation, using the same first-order upwind method, on 
another array of uniform grids, from a coarse one with M a dj = 20 * 2° to a fine 
one with M a dj = 20 * 2 3 . With each of these approximate adjoint solutions, we 
compute an a posteriori estimate of the error in the quantity of interest using 
(18). Then we plot the errors in the quantity of interest against the number of 
grid cells M in Fig. 1. For this simple case, solving the adjoint equation with 
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Figure 1: For the simple kernel function = 1, the errors in the quantity of 
interest for different grids for the primal and adjoint approximations. The same 
first-order method is used for discretizing the primal and adjoint equations. 

the same first-order upwind method on a grid that is at most as fine as the grid 
for the primal equation proves adequate. In fact, with M a dj — 20 * 2 3 (one level 
coarser than the finest grid for solving the primal equation), the a posteriori 
error estimates are indistinguishable from the true errors. 
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However, with a more challenging kernel function, solving the adjoint equa- 
tion with a numerical method of the same order on a grid with similar resolution 
may not be adequate. In such cases, increasing the grid resolution for the ad- 
joint equation can solve this difficulty. A more effective solution strategy is 
to use a higher-order method for the adjoint equation. We demonstrate these 
observations with the choice 

± Mff^,. (24) 

This kernel function emphasizes the accuracy of the solution in a small region 
of radius e surrounding the point (x — L/2,t — T/2) in the space-time compu- 
tational domain. 

We first solve the adjoint equation with the same first-order upwind method 
as used for the primal problem on an array of uniform grids with increasing 
resolutions. The results are presented in Fig. 2. What we see is that on coarse 
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Figure 2: For the kernel function (24), the errors in the quantity of interest for 
different grids for the primal and adjoint approximations. The same first-order 
method is used for discretizing the primal and adjoint equations. 

grids (in this case M ai y = 20 * 2 q with q = 0, 1, 2, 3), the adjoint solutions are 
not accurate enough to produce sensible estimates of the errors in the quan- 
tity of interest. We also see that with increasing resolutions for the adjoint 
approximation, the a posteriori estimates of the errors are converging to the 
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true errors. In fact, at M a dj = 20 * 2 5 , the errors estimates can be regarded as 
decent approximations of the true errors. 

Next, we change our strategy, and solve the adjoint equation using the leap- 
frog method which is second-order accurate both in time and space. The results 
are shown in Fig. 3. The leap-frog method provides a very efficient solution. 




0.5 1 1.5 2 2.5 3 3.5 4 

log2(nCells) 



Figure 3: For the kernel function (24), the errors in the quantity of interest 
for different grids for the primal and adjoint approximations. A higher-order 
method is used for discretizing the adjoint equations compared to that used for 
the primal equations. 

In fact, even on the coarsest grid with just M a dj = 20, the adjoint solution is 
accurate enough to produce error estimates that are indistinguishable from the 
true errors. 

In the Galerkin-type variational approach towards a posteriori error esti- 
mation for finite volume schemes ([4, 9]), it is required that the solution of the 
adjoint problem be sought in a function space larger than the solution space 
for the primal problem, due to the orthogonality between the residual and the 
test function space. This requirement does not apply to our approach, because 
neither of the finite volume schemes for the primal or adjoint problems is based 
on variational formulations. In fact, for the cases (23) and (24) considered in 
this section, both the primal and adjoint problems are solved in piecewise con- 
stant function spaces, so to speak. The more complex case (24) poses some 
challenges because the errors in the adjoint solution starts to deteriorate the a 
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posteriori estimates of the errors in the quantity of interest corresponding to 
this kernel function. Solving the adjoint problem in a broader function space, 
e.g. by implementing a Godunov-type higher-order finite volume scheme ([4, ?]), 
is certainly a possible strategy. This strategy is not explored in this work. In- 
stead, we demonstrate through numerical results that the challenges can also be 
met, to various extents, solely by refining the mesh or implementing a plain high 
order scheme (e.g. from a first-order upwind scheme to a second-order leap-frog 
scheme). The function spaces for the solutions of the adjoint problems remain 
piecewise constant. 

3 Application to adaptive mesh refinement for 
the ID shallow water equations 

In this section, we demonstrate the application of the a priori error estimate, 
derived in Section 2.1, to adaptive mesh refinement for finite volume methods. 
We take the system of one dimensional linearized shallow water equations as an 
example. We also discuss the challenge and the handling of a nonlinear quantity 
of interest. 

3.1 Types of adaptive mesh refinement strategies 

We can identify the following types of adaptive mesh refinement strategies: 

Type 1. Time unvarying adaptive mesh refinement. 

The mesh for the spatial dimensions is adaptively refined, but remains 
fixed for the whole simulation period. The time steps are non-adaptive, 
though they could be nonuniform. A generic depiction of the resulting 
grid is given in Fig. 4. 

Type 2. Non-incremental time varying adaptive mesh refinement. 

The spatial mesh is adaptively refined for each temporal sub-interval of 
the simulation. The time step is non-adaptive, though it could be nonuni- 
form. The mesh refinement is done after each full simulation, i.e., not 
incrementally in time. A generic depiction of the resulting grid is given in 
Fig. 5 

Type 3. Incremental time varying adaptive mesh refinement. 

The spatial mesh is adaptively refined at each sub-interval of the simula- 
tion period, as the simulation proceeds. The time step is non-adaptive, 
though it could be nonuniform. 

Type 4. Incremental time varying adaptive mesh with adaptive time steps. 

In this article, we discuss the Type 1 and 2 strategies and leave the more so- 
phisticated Type 3 and 4 strategies to future work. 

The elegant formula (5) cannot be applied directly for adaptive mesh refine- 
ments because it does not provide information about local errors. In order to 
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Figure 4: Sketch of a Type 1 adaptively refined mesh. 



obtain such information, we need to break down the whole space-time domain 
into slabs, and find a means for calculating the error contribution from each 
slab. We note that each slab can contain one or more of the time steps used to 
for discretization of the partial differential equations. Let {tj}jL be a partition 

of the whole simulation period [0, T], and for each j, let {-KyJ-jJj be a partition 
of the spatial domain (0, 1) during the time period [tj—x,tj]. Then each slab 
Sij is given by 



We note that this partition of the space-time domain can accommodate the 
Type 1 adaptive mesh refinement strategy with N = 1; see Fig. 4. With N > 1, 
it can accommodate the adaptive mesh refinement strategies of Type 2 as well 
as of Type 3 and 4; see Fig. 5. 




By (5), 



Q{u)-Q(u#) = (f-Lu*,u) 

= I (/ - Lu*)udxdt 
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Figure 5: Sketch of a Type 2 adaptively refined mesh. 



where 

Eij = / (/- Lu*)udxdt 

JSij 

denoting the error contribution from each slab Sij. 

Below is a simple adaptive mesh refinement algorithm based on the error 
contribution from each slab: 

Step 0. Specify the total error tolerance TOL. 
Step 1. Calculate E_ij for each i and j . 
Step 2. Calculate Total_error = sum of E_ij 
Step 3. If Total_error >= TOL, then 

For 1 <= j <= N, 1 <= i <= M_j 
If | E_ij | >= TOL/M then 

Refine the grid over the cell K_ij 

End 

End 

Goto Step 1. 
Else if Total.error < TOL, then 
Exit 

End 

We note that this adaptive mesh refinement algorithm usually leads to over 
refinement because it does not account for the error cancellations among grid 
cells. Nevertheless, it is adequate for our purpose to demonstrate the application 
of the a posteriori error estimate (5) to adaptive mesh refinement for finite 
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volume methods. For more discussions on adaptive mesh refinement algorithms, 
see, e.g. [8]. 



3.2 ID shallow water equations 

We now apply the a posteriori error estimate (5) to adaptive mesh refinement 
for the one-dimensional, linearized shallow water system 

— + h x + u x = 0, 

^ o, n 

— + 2h x + u x = 0. 

at 

We remark that, among many things, the shallow water system can be used 
to model tsunami waves. In (25), the system has been non-dimensionalized, 
coefficients involving physical quantities have be replaced by somewhat artificial 
constants, and the mathematically non-essential Coriolis forcing terms have 
been omitted, all for the sake of simplicity. Despite these heavy simplifications, 
the system (25) still retains some interesting physical features, e.g. gravity 
waves (the underlying mechanism of tsunamis [?, ?]), and therefore is adequate 
for the purpose of this section. 

The coefficient matrix of the system (25), given by 




has two eigenvalues 

A+ = 1 + V2 > and A_ = 1 - V2 < 0. 

The system is sometimes referred to as the subcritical mode of the primitive 
equations due to the the opposite signs of the eigenvalues. For a discussion on 
a related problem, see [6]. The eigenvectors of the coefficient matrix form a 
transformation matrix P such that 




We let 



Then, (25) is transformed to 




P- 1 • (26) 



0. (27) 
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We prescribe the upwind boundary conditions for £ and rj given by 



£ = at x = 0, 
rj = at x = 1. 



(28) 



In summary, the primal problem consists of (25) (or (27)) and the boundary 
conditions (28). We let f2 = (0, 1) x (0, T) and u denote the vector of unknowns 
(h, u) and let 

' h t + h x +u x \ 
u t + 2h x + u x 



Lu 



The domain of the operator L is then defined as 

D(L) = {u e L 2 (n), LueL 2 (£l), and u satisfies (28)} , 

where u is related to £ and rj through (26). Therefore, L is an unbounded 
operator in L 2 (f2) with domain D(L). 

We shall next define the domain D{L*) of the adjoint operator. A function u 
belongs to D(L*) if and only if u (Lu, u) is a linear continuous functional 
on D(L) for the norm of L 2 {Q). Let u e D(L) and let u e C°°(fi). Then, 

(Lu,u) = / (ft t + u /i x + h u x )h + (u t + g h x + u u x )u 
Jn 

(■T 



1 „ t 
hh 



T ~ ~ T 
hht + uu 



o 
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where 



P T u. 



(29) 
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Integrating by parts on the space interval we obtain 

fl r i-T 



(Lu, u) = (hh + uu) 
Jo 

J 

Jn 

= / (hh + uu) 
Jo 



t=T 



[ hh t + [ ( 
Jn Jo v 



t=T 



dx + y (\+€£\ x =i - \-rrn\ x = ^ dt- 



/ hh t + uu t + A+££ x + \-r)r\ x 
J a 
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Jo t=T Jo v ' 



hh t + uu t 



\ (hh + uu) dx + (A+££| x= i - \_r]T]\ x= o) dt - 
Jo t=T Jo v J 



~ ~ t d „ 

hh t + uut + u ■ A — — u 
ox 



■ I (hh + uu) 
Jo 

In ( 







dx + 




t=T 





A+££| x= i - A_77J7| X=0 ) 



dt- 



du ^rp d _ , 

dt dx 



For u — > (Lu, u) to be continuous for the L 2 norm, it is necessary that 

( u = at t = T, 

£ = at x = 1, 
rj = at x = 0, 

du iT d „ r 2/^\ 



(30) 



In the above, (£, 77) is the transformation of u defined in (29). Therefore, we 
define the domain for the adjoint operator as 



D(L*) = {u e L 2 (VL) I u satisfies (30)} 



(31) 



For each u e D(L*), 

du T d „ 

The adjoint problem can be stated as follows: 

/or eijery G L 2 (f2), find u £ D(L*) such that 



L*u = 4>. 
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The well-posedness of the adjoint problem can be established by the semigroup 
theory. 

We consider the nonlinear quantity of interest: 

Q{h, u) = I [ hu 2 (32) 

which is the time integral of the kinetic energy. Let h# and u# be the numerical 
approximations to h and u, respectively. Then, the error in the quantity of 
interest can be computed as 

Q{h, u) - Q(/i # , u*) = f \hu 2 - f \h*u* 2 

= ( {h-h*)\u 2 + f (u-u*)\h*(u + u*). 
Jn 2 J n 2 

Therefore, the kernel function is given by 

<t> = fa, \h*(u + u#)y (33) 

We note that the kernel function <j> involves the unknown function u, an issue 
inevitable for nonlinear quantities of interest. In practice, the true values of the 
unknown functions are not available, and therefore have to be approximated. 
To the best of our knowledge, there is yet no rigorous theory guding the choice 
of the approximations. An obvious option, which is also what is usually taken 
in the literature, is to replace the unknowns by their numerical approximations. 
Because the goal here is to evaluate the applicability of our a posteriori error 
estimation approach to adaptive mesh refinement, we content ourselves with this 
option, and leave the more fundamental questions to future endeavor. Thus in 
what follows, we replace u by u# in (33). 

For the system (25), we specify as the initial conditions 

h(x,0) = h (x) = | 1 \ x -^) <e ' 
[ elsewhere, 

and 

u(x,0) = 0. 

The set of initial conditions represents a flow at rest with a bulk of fluid artifi- 
cially raised above the surface, reminiscent of the sudden flow elevation caused 
by an earthquake under the sea. When the flow is allowed to evolve freely, the 
bulk of fluid will split into two smaller wave packets of equal size and travel in 
opposite directions with different speeds. See Fig. 6 for a snapshot of the solu- 
tions h and u. The problem is numerically challenging due to the discontinuous 
data and solutions and to the wave packets moving in different directions at 
different speeds. 
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Figure 6: Snapshot of the solution of (25). The arrows point to the wave travel 
directions; longer arrows indicate greater speeds, but not proportionally. 

We experiment with Type 1 and Type 2 mesh refinement strategies and 
compare the results with that of the uniform mesh refinement strategy. All 
strategies are applied with the same tolerance goal 

TOL = 4.0 x 1CT 4 

for the error. The uniform mesh refinement strategy requires 2,560 cells to 
reach the tolerance goal. The Type 1 strategy, i.e., time-unvarying adaptive 
mesh refinement, requires 1,200 cells to reach the same tolerance goal. The 
final grid is plotted in Fig. 7. We see that the grid is intensely refined over 
the travel ranges of both of the wave packets. The savings in number of cells 
come from the "quiet" regions where there are no wave activities. The refined 
region is shifted towards the right end because the rightward wave is traveling 
at a greater speed and thus has a longer travel range. We then apply the Type 
2 strategy with three temporal sub-intervals over the whole simulation period 
[0, T\, The resulting grid that meets the tolerance is given in Fig. 8 and has 
949 cells for the first sub-interval, 1143 for the second, and 179 for the last, 
with an average of 757 cells for the whole simulation period. The saving in 
the number of cells stems from the fact that when the wave packets are far 
apart, the surrounding regions for them can be refined separately and thus the 
unnecessary refinement in the middle is avoided. 



18 




4 Concluding remarks 

In this article, we present a framework for goal-oriented a posteriori error es- 
timation for finite volume methods. The formulation of the a posteriori error 
estimate is independent of the exact form of the methods and therefore can 
be applied to arbitrary finite volume methods. In this framework, it is not re- 
quired that the adjoint equations be solved in a function space larger than that 
for the forward equation, due to the fact that finite volume methods, generally 
speaking, are not built on Galerkin orthogonality. 

To demonstrate the validity of the a posteriori error estimate, we conduct 
numerical experiments with the one-dimensional linear transport equation. The 
primal equation is solved using the first-order, upwind finite volume method. 
The overall conclusion from these experiments is that the more accurate the ad- 
joint solution is, the more accurate the a posteriori error estimate is. Roughly 
speaking, there are two scenarios one has to deal with in practice. If the quan- 
tity of interest is defined by a simple kernel function, then solving the adjoint 
equation with a numerical method of the same order or accuracy and with a 
mesh with resolution similar as that for the primal equation may prove ade- 
quate. If the quantity of interest is more complex, then the adjoint equation 
may have to be solved on a finer mesh or to be solved by a high-order accurate 
method. We have found that the latter approach is often more efficient. 

An application of the a posteriori error estimate to adaptive mesh refinement 
is also presented. The one-dimensional linearized shallow-water equations are 
taken as an example. The test case involves two wave packets traveling in oppo- 
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Figure 8: The final grid generated by the Type 2 adaptive mesh refinement 
strategy 



site directions with different speeds. This case is numerically very challenging. 
The a posteriori error estimate is found to be effective at guiding various adap- 
tive mesh refinement strategies that lead to grids that are dynamically refined 
according to wave activities. 

We can identify several directions or future work. The impact of the lack 
of accuracy in the adjoint solution on the accuracy of the a posteriori error 
estimate is only experimentally explored in this work. More rigorous analysis 
is warranted for this issue. A related issue is the impact of the substitution of 
the primal solution u by its approximation w& in the kernel function (j>. We did 
not touch upon this issue here, but it is very important and inevitable in cases 
involving nonlinear quantities of interest. 

Another direction future research is to apply the framework of a posteriori 
error estimation, laid out in this article, to regional climate modeling, which is 
the primary motivation for the current work. An emerging approach towards 
climate modeling is to use one global grid over the whole sphere, with local 
refinements over regions of interest, and with smooth transitions between coarse 
and fine regions ([16]). In our opinion, the current mesh refinement strategy used 
in this approach is quite rudimentary in that it only refines over regions that are 
directly of interest. Goal-oriented error estimation is clearly needed to guide a 
more sensible mesh refinement strategy, and thus to control the computational 
errors with regard to the quantities of interest. 
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